% SCRIPT TO TRACK THE POSITION OF LAGRANGIAN PARTICLES IN A 2D CLOSED DOMAIN
% WITH MOVING BOUNDARIES WHICH VARIES AT DIFFERENT SPEEDS (om1, om2)
close all
clear all
Fig_settings
% GRID PARAMETERS
n=160;
N=n;
R=2.0;
mu=1;
epsilon=0.25*2*pi*R/N;
% BOUNDARY PARAMETERS
om=1;
t_fin=2*pi/om;
dt0=2*pi*om/(20);
t_step=round(t_fin/dt0)+1;
t=0;
phi=0;
% FREQUENCY OF THE BOUNDARY
om1=om;
om2=om/5;
% PARTICLE PARAMETER
Nout=4;
xout=R/2*[cos(0) cos(pi/16) cos(pi/8) cos(pi/4)];
yout=R/2*[sin(0) sin(pi/16) sin(pi/8) sin(pi/4)];
xP=zeros([t_step,Nout]);
yP=zeros([t_step,Nout]);
xLagr=xout;
yLagr=yout;
% VISUALIZE BOUNDARY AND PARTICLE POSITION
[xb0,yb0,ub,vb,theta,dtheta]=moving_coord2D_phi(n,R,t,om,phi);
[xb1,yb1,ub,vb,theta,dtheta]=moving_coord2D_phi(n,R,t_fin,om,phi);
% INIZIALIZATION
dt=dt0;
cont=1;
% SELECT AN APPROPRIATE dt, SUCH THAT THE BOUNDARY DOES NOT MOVE TOO FAST
maxom=max([om1 om2]);
[xb,yb,ub,vb,theta,dtheta]=moving_coord2D_phi(n,R,t,maxom,phi);
dt=dt0;
while t<t_fin
    xb_old=xb;
    yb_old=yb;
    t=t+dt;
    [xb,yb,ub,vb,theta,dtheta]=moving_coord2D_phi(n,R,t,maxom,phi);
    distb=sqrt((xb-xb_old).^2+(yb-yb_old).^2);
    if max(distb)>epsilon
       dt=dt/2;
       t=t-dt;
    end
end
% TRACK THE POSITION OF THE PARTICLES IN TIME
dt0=dt;
t=0;
om1=om;
om2=om*2;
while t<t_fin
    step=1;
    [uP,vP,xb,yb,ub,vb]=particle_coord(n,R,t,om,om1,om2,phi,mu,epsilon,xLagr,yLagr,Nout,dt);
    while  (step == 1)        
        %utemp=uLagr+(1./taup)*(uLagr-uP)*dt;
        %vtemp=vLagr+(1./taup)*(vLagr-vP)*dt;
        xtemp=xLagr+uP*dt;
        ytemp=yLagr+vP*dt;
        for j=1:Nout
            [dist(j), idist(j)]=min(sqrt((xb-xtemp(j)).^2+(yb-ytemp(j)).^2));
        end
        [absmin, iabsmin]=min(dist);
        if (absmin<epsilon)
            dt=dt/2
            step=1;
        else
            dt=dt0;
            step=0;
        end
    end 
    xLagr=xtemp;
    yLagr=ytemp;
    t=t+dt;
    cont=cont+1;
    time_stepping(cont)=dt;
    xP(cont,:)=xLagr(:);
    yP(cont,:)=yLagr(:);

    figure(1)
    plot(xLagr,yLagr,'.','Markersize',3)
    hold on
    plot(xb, yb, 'k-.')
    hold off
    title(['t= ' num2str(t)])
    axis([-4 4 -4 4])
  
    figure(2)
    plot(t,xLagr,'.')
    title(['t= ' num2str(t)])
    hold on
end

figure(3)
plot(xb0, yb0, 'k-.')
hold on
plot(xb1, yb1, '-.')
plot(xout, yout, 'r.')
plot(xP(:,1),yP(:,1),'.')
plot(xP(:,2),yP(:,2),'r.')
plot(xP(:,3),yP(:,3),'m.')
plot(xP(:,4),yP(:,4),'g.')
xlabel('X')
ylabel('Y')
axis([-2*R 2*R -2*R 2*R])
grid on



